1 Requirements

We need theses packages

library(readxl) # for reading xl files
library(ggplot2) # for figures
library(CaliCo) # calibration package
library(RobustCalibration) # another calibration package
library(DiceDesign) # package for numerical design of experiments
library(DiceKriging) # package for computing Gaussian process Emulator
library(tidyr) # manipulate data

We load these home-made functions that are the computer models to calibrate

source("functions.R")

We load some data about ball dropping. The data are used in (Silva et al. 2020) can be found on the github repo.

# the sheet number correspond to the type of ball (see the xl file)
don=read_xlsx("Ball_drops_data.xls",sheet = 5)
names(don)=c("drop","time","Height","Velocity")
don$drop=as.factor(don$drop)

We start with a plot

ggplot(don,aes(x=time,y=Height,col=drop))+geom_point()

2 Simple calibration

We’ll keep the data from the first drop.

don1=don[don$drop==1,]

We consider the model (compute code) that provides the height of the ball as a function of time \(x=t\) and with two parameters to calibrate \(\theta=(g,h_0)\) the gravity constant and the initial height. \[(x,\theta)\mapsto h=f(x,\theta)\] We set the calibration with CaliCo under model 1 according to (Carmassi et al. 2019) terminology:

model1 <- model(balldropg,X=don1$time,Yexp = don1$Height,"model1")
print(model1)
Call:
[1] "model1"

With the function:
function (t, theta) 
{
    g = theta[1]
    h0 = theta[2]
    h = -g * 0.5 * t^2 + h0
    h[h < 0] = 0
    return(h)
}

No surrogate is selected

No discrepancy is added
# we provide some parameters to have a plot
model1 %<% list(theta=c(9.8,46),var=1e-1)
model1$plot(x=don1$time)

We can choose the prior distributions on the two parameters in \(\theta\) and on the observational variance (or measurement error):

opt.prior <- list(c(10,1),c(45,5),c(1e-2,1e-0))
type.prior <- c(rep("gaussian",2),"unif")
pr1 <- prior(type.prior,opt.prior)
plot(pr1$Prior1)

plot(pr1$Prior2)

We can then run the calibration and have a look at the results:

opt.estim = list(Ngibbs=10000,Nmh=10000,thetaInit=c(9.8,45,1e-1),
                 r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,Nchains=1,burnIn=2000)
mdfit1 <- calibrate(model1,pr1,opt.estim)
print(mdfit1)
Call:

With the function:
function (t, theta) 
{
    g = theta[1]
    h0 = theta[2]
    h = -g * 0.5 * t^2 + h0
    h[h < 0] = 0
    return(h)
}
<bytecode: 0x56168fa094d0>

Selected model : model1 

Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "38.08%" "24.03%" "28.6%" 

Acceptation rate of the Metropolis Hastings algorithm:
[1] "32.17%"

Maximum a posteriori:
[1]  9.2204626 46.7244840  0.1101157

Mean a posteriori:
[1]  9.1313193 46.5504446  0.1039487
plot(mdfit1,don1$time)

m = mdfit1$output$out$THETA[-(1:1000),]
hist(m[,2])

hist(m[,1])

hist(m[,3])

By using the MAP for \(\theta\), we can obtain predictions and compute the MSE.

ypred = forecast(mdfit1,don1$time)
plot(ypred,don1$time)

mean((ypred$y.new.pred$y-don1$Height)^2)
[1] 0.09271534

In order to test the prediction and the coverage rate, we can use cross validation

mdfitCV <- calibrate(model1,pr1,
                     opt.estim = list(Ngibbs=1000,
                                      Nmh=5000,
                                      thetaInit=c(9.8,45,1e-1),
                                r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,Nchains=1,burnIn=2000),
                     opt.valid = list(type.valid="loo",
                                      nCV=50))
print(mdfitCV)
Call:

With the function:
function (t, theta) 
{
    g = theta[1]
    h0 = theta[2]
    h = -g * 0.5 * t^2 + h0
    h[h < 0] = 0
    return(h)
}
<bytecode: 0x55a30dead190>

Selected model : model1 

Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "63.5%" "17.3%" "20.5%"

Acceptation rate of the Metropolis Hastings algorithm:
[1] "59.76%"

Maximum a posteriori:
[1]  9.2636016 46.7784532  0.1069267

Mean a posteriori:
[1]  9.1392385 46.5589027  0.1019072


Cross validation:
 Method: loo 
  Predicted     Real       Error
1  43.13085 43.13187 0.001017761
2  41.67053 41.89956 0.229031427
3  31.15424 31.54828 0.394039020
4  23.69923 23.83016 0.130929637
5  19.42488 19.44317 0.018283849
6  39.59713 39.91591 0.318785698

RMSE: [1] 0.506053

Cover rate: 
[1] "100%"
mdfitCV$ResultsCV
   Predicted     Real        Error
1   43.13085 43.13187 0.0010177610
2   41.67053 41.89956 0.2290314272
3   31.15424 31.54828 0.3940390202
4   23.69923 23.83016 0.1309296366
5   19.42488 19.44317 0.0182838492
6   39.59713 39.91591 0.3187856982
7   46.55498 46.10805 0.4469310920
8   27.59080 27.89219 0.3013968836
9   16.38248 16.11230 0.2701750090
10  46.52468 46.13628 0.3883948585
11  11.41430 10.83166 0.5826354628
12  32.23738 32.52318 0.2857968186
13  42.26879 42.44439 0.1756022810
14  44.10040 43.89717 0.2032301222
15  34.29827 34.63545 0.3371832491
16  20.85504 20.98674 0.1317023705
17  42.85000 42.93345 0.0834527333
18  17.92037 17.73712 0.1832528099
19  46.39718 46.02302 0.3741640302
20  11.47859 10.83166 0.6469234430
21  44.53033 44.32234 0.2079927645
22  19.42456 19.44317 0.0186116690
23  44.52353 44.32234 0.2011909910
24  33.30193 33.66056 0.3586250895
25  32.23326 32.52318 0.2899216900
26  26.34542 26.51110 0.1656788965
27  42.28382 42.44439 0.1605746881
28  38.82019 39.10350 0.2833092276
29  45.52874 45.28606 0.2426833398
30  42.86560 42.93345 0.0678502795
31  42.85940 42.93345 0.0740591361
32  46.49433 46.05136 0.4429686972
33  27.59683 27.89219 0.2953659406
34  27.59944 27.89219 0.2927532708
35  41.65894 41.89956 0.2406195746
36  38.82299 39.10350 0.2805170876
37  19.43654 19.44317 0.0066320735
38  46.54895 46.10805 0.4408978018
39  36.22515 36.58504 0.3598898127
40  46.56236 46.13628 0.4260795823
41  19.44277 19.44317 0.0003965228
42  13.11581 12.61896 0.4968491463
43  43.63901 43.55704 0.0819716953
44  27.61064 27.89219 0.2815574037
45  32.22386 32.52318 0.2993148857
46  20.91339 20.98674 0.0733543839
47  36.23090 36.58504 0.3541414523
48  44.87945 44.60579 0.2736662923
49  44.89882 44.60579 0.2930376660
50  45.57710 45.28606 0.2910443663

3 Calibration with discrepancy

We propose then to add a discrepancy term in the model which corresponds to model 3 according to (Carmassi et al. 2019)’s terminology.

model3 <- model(balldropg,don1$time,don1$Height,"model3",
                opt.disc = list(kernel.type="gauss"))
model3 %<% list(theta=c(9.8,46),thetaD=c(1e-1,1),var=1e-3)
plot(model3,don1$time)

type.prior <- c(rep("gaussian",2),"unif","unif","unif")
opt.prior <- list(c(10,1),c(40,5),c(1e-2,1e-0),c(0,2),c(1e-2,1e-0))
pr2 <- prior(type.prior,opt.prior)

opt.estim2=list(Ngibbs=1e4,Nmh=1e4,thetaInit=c(10,46,1e-1,0.5,1e-1),
                r=c(0.1,0.1),sig=(diag(c(1,1,1e-1,.5,1e-1)*.2))^2,Nchains=1,burnIn=2000)
mdfit3 <- calibrate(model3,pr2,opt.estim2)
plot(mdfit3,don1$time)
[1] "The computational time might be long to get the output plot"

print(mdfit3)
Call:

With the function:
function (t, theta) 
{
    g = theta[1]
    h0 = theta[2]
    h = -g * 0.5 * t^2 + h0
    h[h < 0] = 0
    return(h)
}
<bytecode: 0x56157fac3c98>

Selected model : model3 

Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "65.66%" "45.17%" "51.22%" "53.36%" "54.74%"

Acceptation rate of the Metropolis Hastings algorithm:
[1] "27.8%"

Maximum a posteriori:
[1]  9.5662782 47.3139099  0.1192344  1.3893981  0.1057993

Mean a posteriori:
[1]  9.22878566 46.43685325  0.10479450  0.74960623  0.07172582
m = mdfit3$output$out$THETA[-(1:1000),] 
hist(m[,1])

hist(m[,2])

hist(m[,3])

hist(m[,4])

hist(m[,5])

By using the MAP for \(\theta\), we can obtain predictions and compute the MSE.

ypred = forecast(mdfit3,don1$time)
plot(ypred,don1$time)

mean((ypred$y.new.pred$y-don1$Height)^2)
[1] 0.208324

4 Time-consuming code

We assume here that the code modeling the ball drop is time-consuming. We make recourse to a Gaussian Process Emulator. We first consider that there is no discrepancy. Hence, we are in the framework of model 2.

# bounds on theta for the GP emulator
binf = c(8,40)
bsup = c(11,50)
# design of experiments and limited number of evaluations to the computer code
ndesign=10
DOE <- maximinSA_LHS(lhsDesign(ndesign,3)$design)$design
DOE <- unscale(DOE,c(min(don1$time),binf),c(max(don1$time),bsup))
Ysim <- balldropg(DOE[1,1],DOE[1,2:3])
for (i in 2:ndesign){Ysim <- c(Ysim,balldropg(DOE[i,1],DOE[i,2:3]))}

We provide the design of numerical experiments with the associated evaluations:

model2 <- model(code=NULL,don1$time,don1$Height,"model2",
                opt.gp = list(type="matern5_2", DOE=DOE),
                opt.sim = list(Ysim=Ysim,DOEsim=DOE))


ParamList <- list(theta=c(9.8,46),var=1e-1)
model2 %<% ParamList

and we plot:

plot(model2,don1$time,CI="all")

plot(model2,don1$time,CI="GP")

opt.estim = list(Ngibbs=1e4,Nmh=1e4,thetaInit=c(9.8,45,1e-1),
                 r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,Nchains=1,burnIn=2000)
mdfit2 <- calibrate(model2,pr1,opt.estim)
mdfit2
Call:

With the function:
NULL

Selected model : model2 

Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "48.48%" "25.08%" "28.01%"

Acceptation rate of the Metropolis Hastings algorithm:
[1] "24.13%"

Maximum a posteriori:
[1] 10.15966591 52.44318131  0.02237512

Mean a posteriori:
[1]  9.97885847 47.97409824  0.02184157
plot(mdfit2,don1$time,graph="result") 

m = mdfit2$output$out$THETA[-(1:1000),] 
hist(m[,1])

hist(m[,2])

hist(m[,3])

4.1 Sequential design of numerical experiments

We try to make additional evaluations of the model with the purpose to improve the calibration.

model2 <- model(balldropg,don1$time,don1$Height,"model2",
                opt.gp = list(type="matern5_2",DOE=NULL),
                opt.emul = list(p=2,n.emul=20,binf=binf,bsup=bsup))


newModel2 <- sequentialDesign(model2,pr1,
                              opt.estim = list(Ngibbs=100,
                                               Nmh=600,
                                               thetaInit=c(9.8,45,1e-1),
                                               r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,
                                               Nchains=1,
                                               burnIn=200),
                              k=10)

We can print the design with added points and make a plot

newModel2$doe.new
          doe.X                   
      2.2437746  8.328901 46.79794
      1.7921889 10.703121 41.68269
      1.3998149  8.741023 42.65614
      1.0054142 10.681355 45.91142
      0.2023622  8.928212 44.82056
      2.1777346 10.228825 48.02350
      0.7633450 10.035511 40.85691
      2.6758861  9.622211 43.49616
      0.4874267  8.293524 49.29969
      1.2030784  9.237339 47.59909
y.new 0.0000000  9.585993 47.48981
y.new 2.7716651  8.771090 46.12897
y.new 0.6683333  8.913166 46.53728
y.new 1.9033322  8.926643 46.61900
y.new 2.4366653  9.025496 46.42469
y.new 0.0000000  8.939712 46.07890
y.new 2.7716651  9.072577 46.50756
y.new 0.2683333  9.072448 46.62499
y.new 1.3683325  9.012115 46.27818
y.new 1.0349994  8.960061 46.26790
y.new 2.0366655  9.119028 46.61820
p <- plot(newModel2,don1$time,graph=FALSE)
p$doe[[1]]

And then rerun the calibration with the new design.

model2new = model(code=NULL,don1$time,don1$Height,"model2",
                  opt.gp = list(type="matern5_2", DOE=newModel2$doe.new),
                  opt.sim = list(Ysim=as.vector(newModel2$GP.new@y),DOEsim=newModel2$doe.new))

mdfit2seq <- calibrate(model2new,pr1,opt.estim)
mdfit2seq
plot(mdfit2seq,don1$time,graph="result") 

m = mdfit2seq$output$out$THETA[-(1:1000),] 
hist(m[,1])

hist(m[,2])

hist(m[,3])

5 Time-consuming code and discrepancy

We are considering model 4. But the calibration is really long…

model4 <- model(code=NULL,don1$time,don1$Height,"model4",
                opt.gp = list(type="matern5_2", DOE=NULL),
                opt.sim = list(Ysim=Ysim,DOEsim=DOE),
                opt.disc = list(kernel.type="gauss"))

optimisation start
------------------
* estimation method   : MLE 
* optimisation method : BFGS 
* analytical gradient : used
* trend model : ~1
* covariance model : 
  - type :  matern5_2 
  - nugget : NO
  - parameters lower bounds :  1e-10 1e-10 1e-10 
  - parameters upper bounds :  4.975254 5.216081 19.53854 
  - best initial criterion value(s) :  -30.49232 

N = 3, M = 5 machine precision = 2.22045e-16
At X0, 0 variables are exactly at the bounds
At iterate     0  f=       30.492  |proj g|=       1.2171
At iterate     1  f =       30.234  |proj g|=       0.10773
At iterate     2  f =       30.197  |proj g|=       0.10669
At iterate     3  f =       30.181  |proj g|=       0.10273
At iterate     4  f =       30.061  |proj g|=       0.49083
At iterate     5  f =       29.969  |proj g|=       0.55312
At iterate     6  f =        29.78  |proj g|=       0.25798
At iterate     7  f =       29.762  |proj g|=      0.072688
At iterate     8  f =        29.76  |proj g|=       0.00743
At iterate     9  f =        29.76  |proj g|=    0.00024236
At iterate    10  f =        29.76  |proj g|=    8.4076e-07

iterations 10
function evaluations 12
segments explored during Cauchy searches 12
BFGS updates skipped 0
active bounds at final generalized Cauchy point 2
norm of the final projected gradient 8.40759e-07
final function value 29.7603

F = 29.7603
final  value 29.760286 
converged
[1] "The surrogate has been set up, you can now use the function"
model4 %<% list(theta=c(9.8,46),thetaD=c(1e-1,0.2),var=1e-1)
plot(model4,don1$time,CI="err")

plot(model4,don1$time,CI="GP")

plot(model4,don1$time,CI="all")

mdfit4 <- calibrate(model4,pr2,opt.estim2)
print(mdfit4)
Call:

With the function:
NULL

Selected model : model4 

Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "71.63%" "58.71%" "65.06%" "67.5%"  "68.54%"

Acceptation rate of the Metropolis Hastings algorithm:
[1] "30.92%"

Maximum a posteriori:
[1]  9.3923176 47.7115764  0.1013969  0.6533405  0.1018994

Mean a posteriori:
[1]  9.10329317 46.46869684  0.08681540  0.53280874  0.07886343
plot(mdfit4,don1$time)
[1] "The computational time might be long to get the output plot"

m = mdfit4$output$out$THETA[-(1:1000),] 
hist(m[,1])

hist(m[,2])

hist(m[,3])

6 Using Scaled GP for the discrepancy

We use the package RobustCalibration that implements the method proposed in (Gu and Wang 2017)

First, we consider a standard GP for the discrepancy function.

model_gasp=rcalibration(design=don1$time, observations=don1$Height, p_theta=2,simul_type=1,
                         math_model=balldropg,theta_range=matrix(c(8,11,40,50),2,2,byrow = TRUE)
                         ,S=10000,S_0=2000,discrepancy_type='GaSP')
2500  of  10000  posterior samples are drawn 
5000  of  10000  posterior samples are drawn 
7500  of  10000  posterior samples are drawn 
10000  of  10000  posterior samples are drawn 
Done with posterior sampling 
3759  of  10000  proposed posterior samples of calibration parameters are accepted 
5128  of  10000  proposed posterior samples of range and nugget parameters are accepted 
m_gasp_pred=predict.rcalibration(model_gasp,don1$time,math_model=balldropg,interval_est=c(0.025,0.975),interval_data = TRUE)
25 percent is completed 
50 percent is completed 
75 percent is completed 
100 percent is completed 
dfgasp = data.frame(time=don1$time,pred=m_gasp_pred@mean,lowIC=m_gasp_pred@interval[,1],
                    upIC=m_gasp_pred@interval[,2],obs=don1$Height,
                    predmathmod=m_gasp_pred@math_model_mean)
dfgasp2 = pivot_longer(dfgasp,cols=2:6,names_to = "type",values_to = "Height")
ggplot()+
  geom_line(data = dfgasp2,aes(x=time,y=Height,col=type,linetype=type))+geom_point(data=don1,aes(x=time,y=Height))

We also compute the MSEs for the prediction from the mathematical model or from the mathematical model corrected by the discrepancy function.

## MSE if the math model and discrepancy are used for prediction
mean((don1$Height-m_gasp_pred@mean)^2)
[1] 0.0007628402
# MSE if the math model is used for prediction
mean((don1$Height-m_gasp_pred@math_model_mean)^2)
[1] 0.1791588

Second, we consider a scaled GP for the discrepancy function

model_sgasp=rcalibration(design=don1$time, observations=don1$Height, p_theta=2,simul_type=1,
                         math_model=balldropg,theta_range=matrix(c(8,11,40,50),2,2,byrow = TRUE)
                         ,S=10000,S_0=2000,discrepancy_type='S-GaSP')
2500  of  10000  posterior samples are drawn 
5000  of  10000  posterior samples are drawn 
7500  of  10000  posterior samples are drawn 
10000  of  10000  posterior samples are drawn 
Done with posterior sampling 
2014  of  10000  proposed posterior samples of calibration parameters are accepted 
5611  of  10000  proposed posterior samples of range and nugget parameters are accepted 
m_sgasp_pred=predict.rcalibration(model_sgasp,don1$time,math_model=balldropg,interval_est=c(0.025,0.975),interval_data = TRUE)
25 percent is completed 
50 percent is completed 
75 percent is completed 
100 percent is completed 
dfsgasp = data.frame(time=don1$time,pred=m_sgasp_pred@mean,lowIC=m_sgasp_pred@interval[,1],
                    upIC=m_sgasp_pred@interval[,2],obs=don1$Height,
                    predmathmod=m_sgasp_pred@math_model_mean)
dfsgasp2 = pivot_longer(dfsgasp,cols=2:6,names_to = "type",values_to = "Height")
ggplot()+
  geom_line(data = dfsgasp2,aes(x=time,y=Height,col=type,linetype=type))+geom_point(data=don1,aes(x=time,y=Height))

We also compute the MSEs for the prediction from the mathematical model or from the mathematical model corrected by the discrepancy function.

## MSE if the math model and discrepancy are used for prediction
mean((don1$Height-m_sgasp_pred@mean)^2)
[1] 0.0008628687
# MSE if the math model is used for prediction
mean((don1$Height-m_sgasp_pred@math_model_mean)^2)
[1] 0.09276917

7 To do

  • Calibrate with the data from different balls.
  • Calibrate the function balldropgfrot that considers an additional fluid drag force.

References

Carmassi, Mathieu, Pierre Barbillon, Matthieu Chiodetti, Merlin Keller, and Eric Parent. 2019. “Bayesian Calibration of a Numerical Code for Prediction.” Journal de La Société Française de Statistique 160 (1): 1–30. http://www.numdam.org/item/JSFS_2019__160_1_1_0/.
Gu, Mengyang, and Long Wang. 2017. “An Improved Approach to Bayesian Computer Model Calibration and Prediction.” arXiv Preprint arXiv:1707.08215.
Silva, Brian M de, David M Higdon, Steven L Brunton, and J Nathan Kutz. 2020. “Discovery of Physics from Data: Universal Laws and Discrepancies.” Frontiers in Artificial Intelligence 3: 25.